clear all; clc; close all;
save_res=1;
 
alpha_s=deg2rad(60);
beta_s=deg2rad(0);
 
beta_j=deg2rad(0);
 
alpha=deg2rad(-180:2:180);
beta=deg2rad(-90:3:90);
 
Fa=ones(length(beta),1)*(1+cos(alpha -pi/2)).^2;
[alpha_mesh,beta_mesh]=meshgrid(alpha,beta);
[xa,ya,za]=sph2cart(beta_mesh,alpha_mesh,Fa);
 
figure(1);
surf(xa,ya,za);
xlabel('X'); ylabel('Y'); zlabel('Z');
title('Radiation pattern for one AM');
axis equal
 
F=nan(length(beta), length(alpha));
 
q_dB=3;
q=10^(q_dB/10);
 
figure(2)
pos=get(gcf,'Position'); pos(3)=800; set(gcf,'Position',pos);
 
for alpha_j=deg2rad(30:5:95);
    
    C=H(alpha_j, beta_j);
    Dn=q*C*C'+eye(2);
    Hf=H(alpha_s, beta_s);
    beta_w=Dn\Hf/(Hf'*(Dn\Hf));
    for a=1:length(alpha);
        for b=1:length(beta);
            U=beta_w'*H(alpha(a),beta(b));
            F(b,a)=abs(U)^2;
        end
    end
    
    F=Fa.*F;
    
    [x,y,z]=sph2cart(beta_mesh,alpha_mesh,F);
    b0=ceil(length(beta)/2);
    Fb0=F(b0,:);
    
    subplot(1,2,1)
    surf(x,y,z);
    xlabel('X'); ylabel('Y'); zlabel('Z');
    axis equal
    
    subplot(1,2,2)
    polar(alpha,Fb0);
    hold on
    polar([alpha_s alpha_s],[0 max(Fb0)],'g');
    polar([alpha_j alpha_j],[0 max(Fb0)],'r');
    hold off
    drawnow
    if save_res
        s=sprintf('pic/DN_alpha_j_%03.0f.png', round(rad2deg(alpha_j)));
        saveas(gcf,s,'png');
        fprintf('Figure is saved at %s\n',s)
    end
end
